options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(1,4))
      drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
      drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
      drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
      drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}

non linear regression with one explanatory variable

fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

ex8-0.stan

single linear regression

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,20)
y=rnorm(n,x*2+5,5)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)
## Initial log joint probability = -5186.51 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       19      -43.7983     0.0117007    0.00153088           1           1       40    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__     -43.80
##    b0         3.31
##    b1         2.13
##    s          5.42
##    m0[1]     11.66
##    m0[2]     22.90
##    m0[3]     42.49
##    m0[4]     11.58
##    m0[5]     19.80
##    m0[6]     39.50
##    m0[7]     26.65
##    m0[8]     13.29
##    m0[9]     21.96
##    m0[10]     9.68
##    m0[11]    24.23
##    m0[12]    31.88
##    m0[13]    26.46
##    m0[14]    40.10
##    m0[15]    33.99
##    m0[16]     4.49
##    m0[17]    40.67
##    m0[18]    19.25
##    m0[19]     6.71
##    m0[20]    18.32
##    y0[1]     14.41
##    y0[2]     14.42
##    y0[3]     44.48
##    y0[4]      4.80
##    y0[5]     20.06
##    y0[6]     41.77
##    y0[7]     31.97
##    y0[8]      3.64
##    y0[9]     17.72
##    y0[10]    13.58
##    y0[11]    21.38
##    y0[12]    32.26
##    y0[13]    21.65
##    y0[14]    35.87
##    y0[15]    42.02
##    y0[16]    -7.79
##    y0[17]    40.27
##    y0[18]    17.12
##    y0[19]     3.28
##    y0[20]    19.37
##    m1[1]      4.49
##    m1[2]      8.72
##    m1[3]     12.94
##    m1[4]     17.16
##    m1[5]     21.38
##    m1[6]     25.60
##    m1[7]     29.83
##    m1[8]     34.05
##    m1[9]     38.27
##    m1[10]    42.49
##    y1[1]     -4.85
##    y1[2]      6.70
##    y1[3]     11.99
##    y1[4]     17.65
##    y1[5]     29.19
##    y1[6]     24.13
##    y1[7]     20.32
##    y1[8]     39.09
##    y1[9]     39.25
##    y1[10]    45.54
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -43.71 -43.36 1.36 1.07 -46.44 -42.27 1.01      563      555
##    b0       3.30   3.18 2.74 2.66  -0.73   7.98 1.02      283      325
##    b1       2.14   2.14 0.25 0.24   1.71   2.54 1.01      367      539
##    s        6.18   6.04 1.15 1.03   4.63   8.35 1.00     1278     1122
##    m0[1]   11.66  11.54 1.96 1.86   8.75  14.93 1.02      299      431
##    m0[2]   22.91  22.89 1.39 1.28  20.61  25.26 1.01      788     1047
##    m0[3]   42.53  42.57 2.68 2.57  38.08  46.80 1.00      859     1286
##    m0[4]   11.58  11.47 1.96 1.88   8.66  14.87 1.02      299      431
##    m0[5]   19.81  19.81 1.45 1.34  17.47  22.28 1.01      496      654
##    m0[6]   39.54  39.60 2.39 2.29  35.60  43.32 1.00     1052     1234
##    m0[7]   26.68  26.67 1.45 1.39  24.26  29.09 1.00     1765     1321
##    m0[8]   13.30  13.19 1.82 1.72  10.53  16.40 1.02      310      435
##    m0[9]   21.97  21.95 1.40 1.31  19.70  24.32 1.01      669      960
##    m0[10]   9.68   9.57 2.13 2.01   6.50  13.29 1.02      292      401
##    m0[11]  24.25  24.26 1.40 1.30  21.90  26.55 1.01     1065     1203
##    m0[12]  31.91  31.94 1.73 1.66  29.01  34.81 1.00     2098     1391
##    m0[13]  26.48  26.47 1.45 1.39  24.08  28.87 1.00     1721     1338
##    m0[14]  40.14  40.19 2.44 2.35  36.07  44.02 1.00     1000     1250
##    m0[15]  34.03  34.06 1.89 1.82  30.90  37.18 1.00     1875     1358
##    m0[16]   4.48   4.37 2.62 2.55   0.61   8.94 1.02      284      325
##    m0[17]  40.71  40.75 2.50 2.38  36.57  44.69 1.00      956     1284
##    m0[18]  19.26  19.26 1.47 1.35  16.89  21.75 1.01      467      627
##    m0[19]   6.71   6.60 2.40 2.32   3.15  10.75 1.02      286      357
##    m0[20]  18.33  18.33 1.51 1.36  15.93  20.92 1.02      426      568
##    y0[1]   11.44  11.42 6.54 6.34   0.90  22.24 1.00     1574     1696
##    y0[2]   22.75  22.71 6.39 6.28  11.82  32.46 1.00     2067     1917
##    y0[3]   42.60  42.66 6.70 6.63  31.83  53.94 1.00     1990     1839
##    y0[4]   11.61  11.56 6.54 6.51   0.57  22.28 1.00     1250     1616
##    y0[5]   20.05  20.00 6.53 6.29   9.09  30.64 1.00     2005     1981
##    y0[6]   39.51  39.62 6.80 6.64  28.57  49.98 1.00     1896     1785
##    y0[7]   26.71  26.83 6.48 6.09  16.03  37.37 1.00     1811     2013
##    y0[8]   13.50  13.34 6.42 6.22   3.32  23.64 1.00     1809     1724
##    y0[9]   22.14  22.11 6.42 6.20  11.51  32.87 1.00     1789     1930
##    y0[10]   9.97   9.97 6.65 6.26  -0.64  21.04 1.00     1286     1727
##    y0[11]  24.55  24.54 6.40 5.99  13.83  35.06 1.00     2118     1735
##    y0[12]  32.04  32.07 6.70 6.52  21.19  43.09 1.00     2099     2025
##    y0[13]  26.49  26.50 6.55 6.23  15.52  37.34 1.00     2088     1922
##    y0[14]  40.34  40.40 6.93 6.54  28.95  51.26 1.00     1701     1914
##    y0[15]  33.87  33.90 6.69 6.44  22.78  44.80 1.00     1893     1902
##    y0[16]   4.56   4.43 6.96 6.54  -6.94  16.39 1.00     1032     1329
##    y0[17]  40.75  40.62 6.98 6.55  29.26  52.18 1.00     1916     1681
##    y0[18]  19.16  19.12 6.32 6.27   8.78  29.44 1.00     1818     1744
##    y0[19]   6.78   6.76 6.47 6.36  -3.50  17.41 1.00     1267     1644
##    y0[20]  18.30  18.19 6.68 6.53   7.77  29.45 1.00     1517     1945
##    m1[1]    4.48   4.37 2.62 2.55   0.61   8.94 1.02      284      325
##    m1[2]    8.71   8.60 2.21 2.09   5.42  12.43 1.02      289      384
##    m1[3]   12.94  12.83 1.85 1.75  10.16  16.05 1.02      307      435
##    m1[4]   17.17  17.14 1.57 1.43  14.68  19.86 1.02      383      408
##    m1[5]   21.40  21.38 1.41 1.32  19.10  23.77 1.01      613      863
##    m1[6]   25.62  25.60 1.42 1.34  23.22  27.97 1.00     1490     1275
##    m1[7]   29.85  29.88 1.60 1.51  27.15  32.50 1.00     2164     1406
##    m1[8]   34.08  34.11 1.90 1.82  30.94  37.25 1.00     1866     1356
##    m1[9]   38.31  38.36 2.27 2.19  34.60  41.95 1.00     1179     1285
##    m1[10]  42.53  42.57 2.68 2.57  38.08  46.80 1.00      859     1286
##    y1[1]    4.58   4.59 6.94 6.61  -6.93  16.18 1.00     1101     1436
##    y1[2]    8.54   8.59 6.68 6.45  -2.33  19.49 1.00     1457     1701
##    y1[3]   12.67  12.50 6.39 5.90   2.17  23.60 1.00     1350     1587
##    y1[4]   16.93  17.03 6.59 6.43   6.17  27.76 1.00     1645     1790
##    y1[5]   21.53  21.58 6.29 6.04  11.10  31.80 1.00     1831     1918
##    y1[6]   25.83  25.86 6.39 6.03  15.49  35.93 1.00     2147     1991
##    y1[7]   29.65  29.71 6.60 6.30  19.23  40.12 1.00     2000     1817
##    y1[8]   34.07  33.95 6.53 6.19  23.09  44.91 1.00     1850     1825
##    y1[9]   38.08  37.96 6.79 6.35  27.22  49.30 1.00     1672     1851
##    y1[10]  42.51  42.50 6.86 6.63  31.28  53.59 1.00     1831     1649

quadratic regression

y=b0+b2(x-b1)**2

ex8-1.stan

quadratic regression

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real b2;
  real<lower=0> s;
}
model{
  y~normal(b0+b2*(x-b1)^2,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b2*(x[i]-b1)^2;
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b2*(x1[i]-b1)^2;
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,0.5*(x-4)**2+5,1)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-1.stan')
fn(mdl,data)
## Initial log joint probability = -89.708 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       28      -4.82174   0.000298523   7.80957e-06           1           1       39    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__      -4.82
##    b0         4.63
##    b1         4.04
##    b2         0.56
##    s          0.77
##    m0[1]      5.18
##    m0[2]      8.11
##    m0[3]      5.53
##    m0[4]      8.86
##    m0[5]      4.63
##    m0[6]      9.80
##    m0[7]      5.15
##    m0[8]     12.27
##    m0[9]      7.66
##    m0[10]    15.36
##    m0[11]     6.99
##    m0[12]    12.21
##    m0[13]     4.88
##    m0[14]     4.95
##    m0[15]     8.12
##    m0[16]     7.16
##    m0[17]     4.94
##    m0[18]     7.08
##    m0[19]    10.05
##    m0[20]    11.07
##    y0[1]      6.13
##    y0[2]      7.25
##    y0[3]      6.19
##    y0[4]      8.25
##    y0[5]      5.36
##    y0[6]      9.58
##    y0[7]      6.34
##    y0[8]     11.21
##    y0[9]      6.26
##    y0[10]    14.09
##    y0[11]     6.99
##    y0[12]    12.24
##    y0[13]     4.74
##    y0[14]     6.66
##    y0[15]     8.64
##    y0[16]     8.00
##    y0[17]     5.07
##    y0[18]     7.06
##    y0[19]    10.41
##    y0[20]    11.04
##    m1[1]     12.21
##    m1[2]      8.97
##    m1[3]      6.62
##    m1[4]      5.18
##    m1[5]      4.63
##    m1[6]      4.98
##    m1[7]      6.23
##    m1[8]      8.37
##    m1[9]     11.42
##    m1[10]    15.36
##    y1[1]     12.68
##    y1[2]      8.86
##    y1[3]      7.05
##    y1[4]      4.56
##    y1[5]      3.89
##    y1[6]      5.40
##    y1[7]      6.13
##    y1[8]      8.43
##    y1[9]     10.26
##    y1[10]    14.49
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 finished in 0.3 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
## 
##  variable  mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##    lp__   -9.25  -7.12 6.83 1.65 -32.40 -5.49 1.11       24       11
##    b0      4.65   4.64 0.49 0.30   4.08  5.18 1.03      309      289
##    b1      2.76   4.04 4.86 0.08 -15.61  4.18 1.09       29       11
##    b2      0.52   0.55 0.14 0.04   0.01  0.62 1.10       27       12
##    s       1.14   0.90 0.86 0.19   0.68  3.76 1.12       37       11
##    m0[1]   5.40   5.21 0.86 0.32   4.71  7.64 1.10       28       11
##    m0[2]   8.16   8.13 0.47 0.32   7.56  8.78 1.05     1028       24
##    m0[3]   5.73   5.57 0.81 0.31   5.07  7.73 1.11       28       11
##    m0[4]   8.73   8.83 0.59 0.33   7.47  9.35 1.08       34       11
##    m0[5]   4.87   4.67 0.92 0.31   4.12  7.32 1.11       28       11
##    m0[6]   9.59   9.76 0.82 0.37   7.48 10.33 1.09       28       11
##    m0[7]   5.37   5.19 0.87 0.32   4.68  7.63 1.11       28       11
##    m0[8]  12.06  12.24 0.99 0.52   9.90 13.07 1.08       34       11
##    m0[9]   7.74   7.68 0.50 0.31   7.16  8.39 1.05       87       16
##    m0[10] 14.94  15.28 1.69 0.71  10.26 16.43 1.09       29       11
##    m0[11]  7.01   7.00 0.37 0.27   6.50  7.51 1.05      697       58
##    m0[12] 11.81  12.14 1.45 0.51   7.27 12.94 1.10       26       11
##    m0[13]  5.12   4.93 0.91 0.32   4.40  7.51 1.11       28       11
##    m0[14]  5.14   4.98 0.78 0.29   4.47  7.11 1.12       30       11
##    m0[15]  8.04   8.10 0.44 0.30   7.32  8.59 1.06       67       16
##    m0[16]  7.27   7.20 0.55 0.30   6.69  8.07 1.07       45       12
##    m0[17]  5.14   4.98 0.79 0.29   4.47  7.11 1.12       30       11
##    m0[18]  7.20   7.11 0.56 0.31   6.61  8.02 1.07       42       12
##    m0[19]  9.83  10.01 0.88 0.38   7.46 10.61 1.09       27       11
##    m0[20] 10.77  11.02 1.14 0.44   7.37 11.70 1.10       26       11
##    y0[1]   5.45   5.24 1.67 0.98   3.57  7.76 1.04       77       15
##    y0[2]   8.16   8.15 1.41 1.04   6.31  9.99 1.05     2090       61
##    y0[3]   5.70   5.59 1.52 0.99   3.79  7.86 1.05      500       21
##    y0[4]   8.72   8.84 1.65 1.03   6.61 10.58 1.06      972       26
##    y0[5]   4.92   4.68 1.78 1.05   2.97  7.21 1.06       64       13
##    y0[6]   9.58   9.73 1.64 1.03   7.54 11.49 1.05      172       19
##    y0[7]   5.38   5.22 1.76 0.98   3.44  7.69 1.06      207       15
##    y0[8]  12.08  12.24 1.69 1.07  10.01 13.94 1.05      459       26
##    y0[9]   7.72   7.67 1.37 0.97   5.97  9.51 1.04     1888      102
##    y0[10] 14.90  15.23 2.21 1.20  11.71 17.18 1.07       40       12
##    y0[11]  7.05   7.00 1.45 0.97   5.24  8.99 1.04     1921       41
##    y0[12] 11.78  12.10 2.14 1.13   9.08 13.77 1.07       58       12
##    y0[13]  5.11   4.96 1.66 1.03   3.16  7.30 1.06      299       16
##    y0[14]  5.13   5.00 1.56 0.99   3.27  7.42 1.06      255       21
##    y0[15]  8.04   8.09 1.49 0.98   6.13  9.87 1.05     1898       52
##    y0[16]  7.26   7.17 1.49 0.97   5.44  9.25 1.04     1758       35
##    y0[17]  5.13   4.99 1.62 0.96   3.27  7.03 1.04      734       21
##    y0[18]  7.16   7.14 1.48 0.99   5.26  9.00 1.05     1968       47
##    y0[19]  9.84   9.99 1.68 1.03   7.74 11.81 1.06      248       18
##    y0[20] 10.80  10.98 1.74 1.09   8.46 12.78 1.06      113       16
##    m1[1]  11.81  12.14 1.45 0.51   7.27 12.94 1.10       26       11
##    m1[2]   8.83   8.93 0.61 0.33   7.45  9.46 1.08       33       11
##    m1[3]   6.67   6.64 0.40 0.27   6.16  7.30 1.06      244       19
##    m1[4]   5.35   5.21 0.72 0.29   4.70  7.04 1.11       31       11
##    m1[5]   4.87   4.67 0.91 0.31   4.12  7.28 1.12       29       11
##    m1[6]   5.21   5.02 0.89 0.32   4.51  7.54 1.11       28       11
##    m1[7]   6.40   6.27 0.69 0.30   5.78  7.79 1.10       31       11
##    m1[8]   8.41   8.38 0.46 0.33   7.81  9.00 1.04     1504       67
##    m1[9]  11.26  11.39 0.82 0.47   9.75 12.15 1.07       41       11
##    m1[10] 14.94  15.28 1.69 0.71  10.26 16.43 1.09       29       11
##    y1[1]  11.81  12.10 2.01 1.08   9.02 13.89 1.05       50       12
##    y1[2]   8.85   8.91 1.52 1.00   6.89 10.71 1.06     1718       27
##    y1[3]   6.70   6.66 1.51 1.00   4.78  8.60 1.05     1346       49
##    y1[4]   5.38   5.24 1.58 1.00   3.53  7.52 1.05      569       20
##    y1[5]   4.86   4.71 1.69 0.99   2.96  6.81 1.05      167       21
##    y1[6]   5.24   5.05 1.67 1.00   3.32  7.46 1.05      115       21
##    y1[7]   6.35   6.26 1.69 1.00   4.53  8.28 1.06      976       24
##    y1[8]   8.46   8.40 1.52 0.99   6.64 10.37 1.05     2098       38
##    y1[9]  11.26  11.39 1.72 1.06   9.15 13.27 1.05      552       27
##    y1[10] 14.91  15.23 2.21 1.14  11.31 17.22 1.07       53       11

both log regression

log y=b0+b1log x x,y>0 y=exp b0 x**b1

ex8-2.stan

both log regression

data{
  int N;
  vector<lower=0>[N] x;
  vector<lower=0>[N] y;
  int N1;
  vector<lower=0>[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~lognormal(b0+b1*log(x),s);
}
generated quantities{
  vector[N] m0;
  vector<lower=0>[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*log(x[i]);
    y0[i]=lognormal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector<lower=0>[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*log(x1[i]);
    y1[i]=lognormal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,10)
y=exp(rnorm(n,log(x)*2+1,1))
grid.arrange(qplot(x,y),
             qplot(log(x),log(y)),
             ncol=2)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -73.2126 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       13      -29.1697    0.00290829   0.000218058           1           1       22    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -29.17
##    b0         1.28
##    b1         1.89
##    s          1.04
##    m0[1]      3.89
##    m0[2]      4.63
##    m0[3]      2.94
##    m0[4]      5.20
##    m0[5]      5.60
##    m0[6]      5.61
##    m0[7]      4.05
##    m0[8]      1.01
##    m0[9]      4.59
##    m0[10]    -2.74
##    m0[11]     2.77
##    m0[12]     3.99
##    m0[13]     5.48
##    m0[14]     4.90
##    m0[15]     5.26
##    m0[16]     4.64
##    m0[17]     5.14
##    m0[18]     5.31
##    m0[19]     5.14
##    m0[20]     5.48
##    y0[1]     39.80
##    y0[2]     18.29
##    y0[3]      9.75
##    y0[4]    217.72
##    y0[5]     63.03
##    y0[6]    509.49
##    y0[7]     85.65
##    y0[8]      5.13
##    y0[9]     53.25
##    y0[10]     0.07
##    y0[11]    42.10
##    y0[12]    18.52
##    y0[13]    87.21
##    y0[14]    36.32
##    y0[15]   124.60
##    y0[16]   111.94
##    y0[17]   135.48
##    y0[18]   465.83
##    y0[19]    53.36
##    y0[20]    17.81
##    m1[1]     -2.74
##    m1[2]      1.64
##    m1[3]      2.85
##    m1[4]      3.58
##    m1[5]      4.11
##    m1[6]      4.52
##    m1[7]      4.86
##    m1[8]      5.14
##    m1[9]      5.39
##    m1[10]     5.61
##    y1[1]      0.01
##    y1[2]      7.92
##    y1[3]     14.91
##    y1[4]    142.32
##    y1[5]    157.01
##    y1[6]    172.96
##    y1[7]    157.79
##    y1[8]     62.92
##    y1[9]    260.96
##    y1[10]   307.46
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "s"    "m0"   "y0"   "m1"   "y1"
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median      sd    mad     q5     q95 rhat ess_bulk ess_tail
##    lp__   -30.75 -30.40    1.32   1.09 -33.32  -29.31 1.00      891      933
##    b0       1.26   1.27    0.46   0.45   0.50    2.00 1.00      488      595
##    b1       1.90   1.90    0.26   0.25   1.49    2.35 1.00      610      738
##    s        1.18   1.15    0.21   0.20   0.89    1.56 1.00     1299     1331
##    m0[1]    3.89   3.89    0.27   0.27   3.44    4.32 1.00      941     1184
##    m0[2]    4.63   4.64    0.28   0.28   4.15    5.09 1.00     1536     1275
##    m0[3]    2.93   2.94    0.31   0.30   2.43    3.43 1.00      597      837
##    m0[4]    5.21   5.22    0.32   0.30   4.67    5.73 1.01     2031     1355
##    m0[5]    5.62   5.62    0.35   0.33   5.03    6.19 1.00     1749     1214
##    m0[6]    5.62   5.63    0.35   0.33   5.04    6.20 1.00     1746     1214
##    m0[7]    4.05   4.05    0.27   0.27   3.59    4.49 1.00     1049     1078
##    m0[8]    0.99   1.00    0.49   0.48   0.17    1.77 1.00      486      595
##    m0[9]    4.59   4.59    0.28   0.27   4.12    5.05 1.00     1497     1226
##    m0[10]  -2.78  -2.75    0.97   0.93  -4.44   -1.22 1.00      513      609
##    m0[11]   2.76   2.76    0.32   0.31   2.23    3.27 1.00      569      764
##    m0[12]   3.99   3.99    0.27   0.27   3.53    4.42 1.00     1006     1126
##    m0[13]   5.49   5.49    0.34   0.32   4.93    6.05 1.00     1853     1250
##    m0[14]   4.91   4.91    0.30   0.29   4.40    5.40 1.01     1861     1317
##    m0[15]   5.27   5.27    0.32   0.31   4.71    5.80 1.00     2008     1355
##    m0[16]   4.65   4.65    0.28   0.28   4.17    5.11 1.00     1550     1388
##    m0[17]   5.15   5.16    0.31   0.30   4.62    5.67 1.01     2026     1375
##    m0[18]   5.32   5.32    0.33   0.31   4.76    5.86 1.00     1978     1356
##    m0[19]   5.15   5.16    0.31   0.30   4.62    5.66 1.01     2027     1375
##    m0[20]   5.49   5.50    0.34   0.32   4.93    6.05 1.00     1849     1250
##    y0[1]  108.47  46.43  345.78  46.27   6.19  350.67 1.00     1979     1958
##    y0[2]  229.84 107.34  454.42 109.12  12.80  819.31 1.00     1982     1877
##    y0[3]   42.76  18.47  113.03  18.91   2.38  140.63 1.00     1991     1932
##    y0[4]  429.60 184.63  987.56 188.28  25.65 1466.08 1.00     1789     1737
##    y0[5]  603.81 273.66 1149.40 276.96  32.56 2119.94 1.00     1886     1506
##    y0[6]  600.09 269.02 1206.65 272.69  38.28 2078.40 1.00     2196     1873
##    y0[7]  113.56  56.22  196.14  55.03   7.77  406.47 1.00     1832     2014
##    y0[8]    6.98   2.83   20.30   2.98   0.34   22.00 1.00     1319     1670
##    y0[9]  225.14 100.24  513.38  99.04  14.03  758.47 1.00     1882     1820
##    y0[10]   0.25   0.06    1.37   0.08   0.01    0.80 1.00      897     1299
##    y0[11]  34.37  15.62   71.60  15.72   2.06  110.92 1.00     1794     1908
##    y0[12] 113.78  50.97  228.42  52.88   6.44  392.13 1.00     1783     1839
##    y0[13] 631.91 232.53 2658.28 233.96  31.53 1979.47 1.00     2015     1866
##    y0[14] 276.54 134.94  491.43 134.66  17.43  922.66 1.00     1981     1814
##    y0[15] 438.80 189.07 1341.15 191.20  27.30 1389.58 1.00     1683     1786
##    y0[16] 249.70 102.40  706.75 106.67  13.23  823.03 1.00     2027     1907
##    y0[17] 355.59 170.31  632.65 171.38  22.80 1226.86 1.00     2205     1952
##    y0[18] 444.03 197.32  997.61 201.42  26.88 1564.43 1.00     2146     1870
##    y0[19] 346.25 159.91  611.51 162.48  21.89 1276.36 1.00     2270     1974
##    y0[20] 543.73 232.44 1930.13 231.14  33.79 1597.37 1.00     1993     1795
##    m1[1]   -2.78  -2.75    0.97   0.93  -4.44   -1.22 1.00      513      609
##    m1[2]    1.62   1.63    0.42   0.42   0.92    2.31 1.00      493      635
##    m1[3]    2.84   2.85    0.31   0.31   2.33    3.35 1.00      582      832
##    m1[4]    3.58   3.58    0.28   0.27   3.12    4.01 1.00      783     1141
##    m1[5]    4.11   4.11    0.27   0.28   3.65    4.54 1.00     1091     1157
##    m1[6]    4.52   4.53    0.28   0.28   4.05    4.98 1.00     1437     1243
##    m1[7]    4.86   4.87    0.30   0.28   4.36    5.35 1.00     1803     1346
##    m1[8]    5.15   5.16    0.31   0.30   4.62    5.67 1.01     2028     1375
##    m1[9]    5.40   5.41    0.33   0.31   4.84    5.95 1.00     1925     1336
##    m1[10]   5.62   5.63    0.35   0.33   5.04    6.20 1.00     1746     1214
##    y1[1]    0.20   0.06    0.89   0.07   0.00    0.79 1.00      911     1193
##    y1[2]   12.32   5.01   48.24   4.94   0.59   36.37 1.00     1644     1829
##    y1[3]   39.36  17.96   95.58  17.90   2.36  125.22 1.00     1728     1722
##    y1[4]   84.15  35.97  209.74  35.62   4.69  297.45 1.00     1924     1506
##    y1[5]  123.15  60.15  207.89  59.94   7.91  426.91 1.00     1971     1856
##    y1[6]  235.77  91.10 1683.98  91.17  12.89  703.81 1.00     2079     1886
##    y1[7]  287.76 129.24  720.95 129.21  15.75  961.33 1.00     2034     1776
##    y1[8]  385.82 181.41  863.84 186.17  22.73 1348.83 1.00     1966     1779
##    y1[9]  510.92 223.51 1139.13 225.45  31.41 1705.22 1.00     1983     1703
##    y1[10] 647.64 280.10 2096.45 282.79  35.76 2112.73 1.00     1875     1856

y0=mcmc$draws('y0')
smy0=summary(y0)

qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(log(y)-log(smy0$median),xlab='obs.-prd.',main='residual')
density(log(y)-log(smy0$median)) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=exp(ml5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=exp(mu5)),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=exp(m)),xy,col='black')

qplot(log(x),log(y),col=I('red'))+
  geom_line(aes(x=log(x),y=ml5),xy,col='darkgray')+
  geom_line(aes(x=log(x),y=mu5),xy,col='darkgray')+
  geom_line(aes(x=log(x),y=log(yl5)),xy,col='lightgray')+
  geom_line(aes(x=log(x),y=log(yu5)),xy,col='lightgray')+
  geom_line(aes(x=log(x),y=m),xy,col='black')

exponential increasing/decreasing

y=b0* exp b1* x -> y~N(b0* exp(b1*x),s)  
log y=log b0+b1* x  -> y~logN(log b0 +b1*x,s)
x,y>0,b0>0
(x=0,y=b0)  
b1>0 x->Infinity,y->Infinity  
b1<0 x->Infinity,y->+0  
n=20
x=runif(n,0,5)
y=rnorm(n,10*exp(-2*x),0.5)
grid.arrange(qplot(x,y),
             qplot(x,log(y)),
             ncol=2)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

ex8-3-1.stan

y=b0* exp b1* x -> y~N(b0* exp(b1*x),s)

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0> b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0*exp(b1*x),s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0*exp(b1*x[i]);
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0*exp(b1*x1[i]);
    y1[i]=normal_rng(m1[i],s);
  }
}
mdl=cmdstan_model('./ex8-3-1.stan')
fn(mdl,data)
## Initial log joint probability = -29.9215 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       15       7.95156   3.54027e-05    0.00348442      0.6984      0.6984       29    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__       7.95
##    b0        10.14
##    b1        -1.93
##    s          0.41
##    m0[1]      0.15
##    m0[2]      0.16
##    m0[3]      1.95
##    m0[4]      3.33
##    m0[5]      0.03
##    m0[6]      6.49
##    m0[7]      8.86
##    m0[8]      0.08
##    m0[9]      0.03
##    m0[10]     0.17
##    m0[11]     0.01
##    m0[12]     0.00
##    m0[13]     0.02
##    m0[14]     0.00
##    m0[15]     0.01
##    m0[16]     0.22
##    m0[17]     0.44
##    m0[18]     0.01
##    m0[19]     1.50
##    m0[20]     0.00
##    y0[1]     -0.41
##    y0[2]      0.24
##    y0[3]      1.58
##    y0[4]      3.18
##    y0[5]      0.34
##    y0[6]      6.59
##    y0[7]      8.59
##    y0[8]      0.19
##    y0[9]      0.01
##    y0[10]     0.10
##    y0[11]    -0.21
##    y0[12]     0.53
##    y0[13]    -0.04
##    y0[14]    -0.08
##    y0[15]     0.20
##    y0[16]     0.01
##    y0[17]     0.11
##    y0[18]    -0.56
##    y0[19]     2.07
##    y0[20]     0.01
##    m1[1]      8.86
##    m1[2]      3.26
##    m1[3]      1.20
##    m1[4]      0.44
##    m1[5]      0.16
##    m1[6]      0.06
##    m1[7]      0.02
##    m1[8]      0.01
##    m1[9]      0.00
##    m1[10]     0.00
##    y1[1]      9.42
##    y1[2]      4.28
##    y1[3]      0.88
##    y1[4]     -0.16
##    y1[5]      0.57
##    y1[6]     -0.08
##    y1[7]      0.23
##    y1[8]      0.38
##    y1[9]     -0.88
##    y1[10]    -0.17
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##    lp__    9.23   9.63 1.46 1.12  6.28 10.77 1.01      535      778
##    b0     10.19  10.19 0.56 0.55  9.27 11.12 1.01     1367      835
##    b1     -1.96  -1.96 0.19 0.18 -2.29 -1.66 1.00     1502     1079
##    s       0.47   0.45 0.09 0.07  0.35  0.63 1.03      331      393
##    m0[1]   0.15   0.14 0.06 0.05  0.07  0.25 1.00     1605     1134
##    m0[2]   0.16   0.15 0.06 0.05  0.07  0.27 1.00     1606     1134
##    m0[3]   1.92   1.92 0.26 0.24  1.50  2.36 1.00     1847     1212
##    m0[4]   3.28   3.29 0.29 0.26  2.80  3.74 1.00     2010     1240
##    m0[5]   0.03   0.03 0.02 0.01  0.01  0.06 1.00     1572     1138
##    m0[6]   6.47   6.47 0.28 0.26  6.02  6.93 1.00     2202     1452
##    m0[7]   8.88   8.88 0.42 0.42  8.18  9.57 1.01     1497      974
##    m0[8]   0.08   0.08 0.04 0.03  0.03  0.15 1.00     1589     1135
##    m0[9]   0.04   0.03 0.02 0.02  0.01  0.07 1.00     1575     1139
##    m0[10]  0.17   0.16 0.06 0.06  0.08  0.28 1.00     1608     1134
##    m0[11]  0.01   0.01 0.01 0.01  0.00  0.03 1.00     1561     1158
##    m0[12]  0.00   0.00 0.00 0.00  0.00  0.01 1.00     1549     1104
##    m0[13]  0.02   0.02 0.01 0.01  0.01  0.05 1.00     1568     1138
##    m0[14]  0.00   0.00 0.00 0.00  0.00  0.00 1.00     1545     1104
##    m0[15]  0.01   0.01 0.01 0.00  0.00  0.02 1.00     1556     1158
##    m0[16]  0.22   0.21 0.08 0.07  0.11  0.36 1.00     1618     1113
##    m0[17]  0.44   0.43 0.12 0.11  0.26  0.66 1.00     1652     1127
##    m0[18]  0.01   0.00 0.00 0.00  0.00  0.01 1.00     1555     1104
##    m0[19]  1.48   1.47 0.24 0.22  1.10  1.88 1.00     1781     1231
##    m0[20]  0.00   0.00 0.00 0.00  0.00  0.00 1.00     1545     1104
##    y0[1]   0.16   0.15 0.49 0.46 -0.62  0.98 1.00     1955     1857
##    y0[2]   0.18   0.16 0.48 0.47 -0.60  0.95 1.00     2157     1888
##    y0[3]   1.91   1.91 0.55 0.52  1.01  2.78 1.00     1869     1923
##    y0[4]   3.28   3.29 0.55 0.51  2.38  4.19 1.00     1609     1493
##    y0[5]   0.04   0.04 0.47 0.45 -0.74  0.83 1.00     1952     1541
##    y0[6]   6.46   6.45 0.56 0.52  5.53  7.34 1.00     2212     1592
##    y0[7]   8.86   8.87 0.62 0.58  7.84  9.85 1.01     1712     1510
##    y0[8]   0.08   0.09 0.48 0.45 -0.72  0.91 1.00     1888     1615
##    y0[9]   0.02   0.04 0.48 0.45 -0.73  0.78 1.00     2053     1605
##    y0[10]  0.17   0.17 0.48 0.47 -0.60  0.98 1.00     2039     1873
##    y0[11]  0.01   0.01 0.49 0.48 -0.80  0.80 1.00     1891     1739
##    y0[12]  0.01   0.00 0.47 0.45 -0.75  0.79 1.00     1797     1704
##    y0[13]  0.03   0.04 0.48 0.44 -0.75  0.81 1.00     1852     1511
##    y0[14]  0.00   0.00 0.46 0.45 -0.78  0.74 1.00     1913     1935
##    y0[15]  0.01   0.00 0.48 0.46 -0.75  0.81 1.00     2118     1806
##    y0[16]  0.23   0.23 0.49 0.45 -0.57  1.01 1.00     1984     1679
##    y0[17]  0.45   0.44 0.49 0.47 -0.37  1.25 1.00     1915     1908
##    y0[18]  0.02   0.02 0.48 0.45 -0.74  0.80 1.00     2028     1803
##    y0[19]  1.46   1.46 0.53 0.51  0.59  2.34 1.00     1888     1749
##    y0[20] -0.01  -0.02 0.48 0.45 -0.81  0.77 1.00     2123     1853
##    m1[1]   8.88   8.88 0.42 0.42  8.18  9.57 1.01     1497      974
##    m1[2]   3.21   3.22 0.29 0.26  2.73  3.67 1.00     1997     1231
##    m1[3]   1.17   1.17 0.21 0.20  0.84  1.55 1.00     1742     1246
##    m1[4]   0.43   0.42 0.12 0.11  0.25  0.65 1.00     1652     1127
##    m1[5]   0.16   0.15 0.06 0.05  0.08  0.27 1.00     1607     1134
##    m1[6]   0.06   0.06 0.03 0.02  0.02  0.12 1.00     1585     1135
##    m1[7]   0.02   0.02 0.01 0.01  0.01  0.05 1.00     1569     1138
##    m1[8]   0.01   0.01 0.01 0.00  0.00  0.02 1.00     1558     1158
##    m1[9]   0.00   0.00 0.00 0.00  0.00  0.01 1.00     1551     1104
##    m1[10]  0.00   0.00 0.00 0.00  0.00  0.00 1.00     1545     1104
##    y1[1]   8.87   8.87 0.61 0.59  7.87  9.86 1.00     1803     1654
##    y1[2]   3.21   3.21 0.55 0.52  2.28  4.11 1.00     2248     1656
##    y1[3]   1.17   1.17 0.53 0.51  0.31  2.02 1.00     2003     1565
##    y1[4]   0.45   0.44 0.49 0.47 -0.30  1.24 1.00     1970     1784
##    y1[5]   0.17   0.16 0.46 0.46 -0.56  0.94 1.00     1934     1879
##    y1[6]   0.04   0.03 0.48 0.47 -0.72  0.86 1.00     1651     1800
##    y1[7]   0.02   0.01 0.47 0.43 -0.74  0.80 1.00     2175     1817
##    y1[8]   0.00   0.01 0.47 0.45 -0.78  0.76 1.00     2146     1525
##    y1[9]   0.00  -0.02 0.46 0.43 -0.74  0.77 1.00     2057     1933
##    y1[10]  0.00  -0.01 0.48 0.46 -0.76  0.79 1.00     1853     1846

ex8-3-2.stan

log y=log b0+b1* x -> y~logN(log b0 +b1*x,s)

data{
  int N;
  vector<lower=0>[N] x;
  vector<lower=0>[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0> b0;
  real<lower=-10,upper=10> b1;
  real<lower=0> s;
}
model{
  y~lognormal(log(b0)+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=log(b0)+b1*x[i];
    y0[i]=lognormal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=log(b0)+b1*x1[i];
    y1[i]=lognormal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,5)
y=rlnorm(n,log(10)-2*x,0.5)
grid.arrange(qplot(x,y),
             qplot(x,log(y)),
             ncol=2)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-3-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -233.923 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       24      -13.9674   0.000121639   0.000125591      0.9884      0.9884       37    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
mle
##  variable estimate
##    lp__     -13.97
##    b0         6.63
##    b1        -1.91
##    s          0.49
##    m0[1]     -7.55
##    m0[2]     -4.26
##    m0[3]     -6.12
##    m0[4]     -3.96
##    m0[5]     -3.09
##    m0[6]     -6.58
##    m0[7]     -0.27
##    m0[8]     -5.69
##    m0[9]     -6.76
##    m0[10]    -0.91
##    m0[11]    -2.66
##    m0[12]    -6.13
##    m0[13]    -3.52
##    m0[14]    -1.36
##    m0[15]    -7.24
##    m0[16]     1.31
##    m0[17]    -5.29
##    m0[18]    -1.67
##    m0[19]    -4.83
##    m0[20]    -3.66
##    y0[1]      0.00
##    y0[2]      0.01
##    y0[3]      0.00
##    y0[4]      0.03
##    y0[5]      0.07
##    y0[6]      0.00
##    y0[7]      0.25
##    y0[8]      0.00
##    y0[9]      0.00
##    y0[10]     0.44
##    y0[11]     0.02
##    y0[12]     0.00
##    y0[13]     0.03
##    y0[14]     0.16
##    y0[15]     0.00
##    y0[16]     3.13
##    y0[17]     0.01
##    y0[18]     0.15
##    y0[19]     0.00
##    y0[20]     0.03
##    m1[1]      1.31
##    m1[2]      0.33
##    m1[3]     -0.66
##    m1[4]     -1.64
##    m1[5]     -2.63
##    m1[6]     -3.61
##    m1[7]     -4.60
##    m1[8]     -5.58
##    m1[9]     -6.57
##    m1[10]    -7.55
##    y1[1]      4.33
##    y1[2]      4.23
##    y1[3]      0.52
##    y1[4]      0.35
##    y1[5]      0.17
##    y1[6]      0.02
##    y1[7]      0.01
##    y1[8]      0.01
##    y1[9]      0.00
##    y1[10]     0.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "s"    "m0"   "y0"   "m1"   "y1"
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -12.86 -12.51 1.42 1.16 -15.69 -11.35 1.00      612      600
##    b0       7.84   7.28 2.95 2.32   4.39  12.83 1.01      330      452
##    b1      -1.94  -1.94 0.10 0.10  -2.11  -1.78 1.01      378      605
##    s        0.56   0.54 0.11 0.10   0.41   0.76 1.00      685      554
##    m0[1]   -7.59  -7.58 0.23 0.22  -7.97  -7.25 1.00      771      964
##    m0[2]   -4.25  -4.25 0.12 0.12  -4.46  -4.05 1.00     1675     1206
##    m0[3]   -6.13  -6.13 0.17 0.16  -6.42  -5.87 1.00     1327     1178
##    m0[4]   -3.94  -3.94 0.12 0.12  -4.15  -3.75 1.00     1165     1363
##    m0[5]   -3.06  -3.06 0.13 0.12  -3.28  -2.85 1.00      587     1047
##    m0[6]   -6.60  -6.60 0.18 0.18  -6.91  -6.32 1.00     1053     1179
##    m0[7]   -0.19  -0.20 0.23 0.22  -0.56   0.20 1.01      335      395
##    m0[8]   -5.70  -5.70 0.15 0.15  -5.96  -5.46 1.00     1719     1210
##    m0[9]   -6.79  -6.78 0.19 0.19  -7.11  -6.49 1.00      979     1137
##    m0[10]  -0.85  -0.85 0.21 0.20  -1.17  -0.51 1.01      344      466
##    m0[11]  -2.62  -2.62 0.14 0.13  -2.85  -2.39 1.00      469      747
##    m0[12]  -6.15  -6.15 0.17 0.16  -6.44  -5.89 1.00     1312     1178
##    m0[13]  -3.50  -3.50 0.13 0.12  -3.70  -3.30 1.00      811     1254
##    m0[14]  -1.30  -1.30 0.19 0.18  -1.60  -0.99 1.01      354      503
##    m0[15]  -7.27  -7.27 0.21 0.21  -7.63  -6.95 1.00      838     1031
##    m0[16]   1.41   1.40 0.31 0.30   0.94   1.91 1.01      330      396
##    m0[17]  -5.29  -5.29 0.14 0.13  -5.53  -5.07 1.00     1827     1328
##    m0[18]  -1.62  -1.62 0.17 0.16  -1.90  -1.33 1.01      367      520
##    m0[19]  -4.83  -4.82 0.13 0.12  -5.05  -4.62 1.00     1888     1290
##    m0[20]  -3.64  -3.64 0.12 0.12  -3.85  -3.45 1.00      903     1224
##    y0[1]    0.00   0.00 0.00 0.00   0.00   0.00 1.00     1687     1848
##    y0[2]    0.02   0.01 0.01 0.01   0.01   0.04 1.00     1835     1739
##    y0[3]    0.00   0.00 0.00 0.00   0.00   0.01 1.00     2031     1879
##    y0[4]    0.02   0.02 0.02 0.01   0.01   0.05 1.00     1832     1871
##    y0[5]    0.06   0.05 0.04 0.02   0.02   0.12 1.00     2057     1999
##    y0[6]    0.00   0.00 0.00 0.00   0.00   0.00 1.00     1985     1786
##    y0[7]    1.00   0.84 0.71 0.49   0.31   2.19 1.00     1175     1568
##    y0[8]    0.00   0.00 0.00 0.00   0.00   0.01 1.00     2046     1854
##    y0[9]    0.00   0.00 0.00 0.00   0.00   0.00 1.00     1708     1809
##    y0[10]   0.51   0.42 0.40 0.23   0.16   1.17 1.01     1293     1425
##    y0[11]   0.09   0.07 0.06 0.04   0.03   0.19 1.00     1943     1524
##    y0[12]   0.00   0.00 0.00 0.00   0.00   0.01 1.00     1860     1823
##    y0[13]   0.04   0.03 0.02 0.02   0.01   0.08 1.00     1899     1900
##    y0[14]   0.33   0.27 0.24 0.15   0.10   0.73 1.00     1729     1624
##    y0[15]   0.00   0.00 0.00 0.00   0.00   0.00 1.00     1539     1768
##    y0[16]   5.12   3.98 4.64 2.35   1.40  11.88 1.00      860     1257
##    y0[17]   0.01   0.01 0.01 0.00   0.00   0.01 1.00     1941     1746
##    y0[18]   0.24   0.19 0.18 0.11   0.08   0.54 1.00     1814     1778
##    y0[19]   0.01   0.01 0.01 0.00   0.00   0.02 1.00     1889     1579
##    y0[20]   0.03   0.03 0.02 0.01   0.01   0.07 1.00     1959     1778
##    m1[1]    1.41   1.40 0.31 0.30   0.94   1.91 1.01      330      396
##    m1[2]    0.41   0.40 0.26 0.25   0.00   0.84 1.01      332      390
##    m1[3]   -0.59  -0.60 0.22 0.21  -0.93  -0.23 1.01      339      437
##    m1[4]   -1.59  -1.59 0.18 0.17  -1.87  -1.30 1.01      365      520
##    m1[5]   -2.59  -2.59 0.14 0.13  -2.82  -2.36 1.00      464      712
##    m1[6]   -3.59  -3.59 0.13 0.12  -3.80  -3.39 1.00      868     1188
##    m1[7]   -4.59  -4.59 0.13 0.12  -4.81  -4.38 1.00     1869     1238
##    m1[8]   -5.59  -5.59 0.15 0.15  -5.84  -5.35 1.00     1753     1258
##    m1[9]   -6.59  -6.59 0.18 0.18  -6.90  -6.30 1.00     1059     1179
##    m1[10]  -7.59  -7.58 0.23 0.22  -7.97  -7.25 1.00      771      964
##    y1[1]    5.30   4.18 8.12 2.50   1.49  12.22 1.00      916     1557
##    y1[2]    1.87   1.52 1.46 0.90   0.58   4.27 1.00      929     1085
##    y1[3]    0.67   0.56 0.50 0.32   0.20   1.46 1.00     1398     1893
##    y1[4]    0.25   0.21 0.17 0.11   0.08   0.54 1.00     1585     1551
##    y1[5]    0.09   0.08 0.06 0.04   0.03   0.20 1.00     1863     1951
##    y1[6]    0.03   0.03 0.02 0.01   0.01   0.07 1.00     1618     1508
##    y1[7]    0.01   0.01 0.01 0.01   0.00   0.02 1.00     2166     1862
##    y1[8]    0.00   0.00 0.00 0.00   0.00   0.01 1.00     1862     1796
##    y1[9]    0.00   0.00 0.00 0.00   0.00   0.00 1.00     1861     1691
##    y1[10]   0.00   0.00 0.00 0.00   0.00   0.00 1.00     1702     1703

y0=mcmc$draws('y0')
smy0=summary(y0)

qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,log(y),col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=log(yl5)),xy,col='lightgray')+
  geom_line(aes(x=x,y=log(yu5)),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

growth curve

y=b0* (1-exp(-b1* x)) -> y~N(1-exp(-b1*x),s)  
x,y>0, b0,b1>0
(x=0,y=0), (x->Infinity,y->b0)

ex8-4.stan

growth curve

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0,upper=100> b0;
  real<lower=0,upper=10> b1;
  real<lower=0> s;
}
model{
  y~normal(b0*(1-exp(-b1*x)),s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0*(1-exp(-b1*x[i]));
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0*(1-exp(-b1*x1[i]));
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,10*(1-exp(-0.5*x)),1)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-4.stan')
fn(mdl,data)
## Initial log joint probability = -123.997 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       17       -12.224   0.000253145   9.03431e-05           1           1       26    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__     -12.22
##    b0        10.23
##    b1         0.51
##    s          1.12
##    m0[1]      5.93
##    m0[2]      7.06
##    m0[3]      9.87
##    m0[4]      7.04
##    m0[5]      3.32
##    m0[6]      9.02
##    m0[7]      2.24
##    m0[8]      9.44
##    m0[9]      9.88
##    m0[10]     9.19
##    m0[11]     8.77
##    m0[12]    10.05
##    m0[13]     9.97
##    m0[14]     9.99
##    m0[15]     9.28
##    m0[16]     8.72
##    m0[17]     9.78
##    m0[18]     9.81
##    m0[19]     5.90
##    m0[20]     9.96
##    y0[1]      6.95
##    y0[2]      7.06
##    y0[3]      8.62
##    y0[4]      8.13
##    y0[5]      3.03
##    y0[6]      9.95
##    y0[7]      1.84
##    y0[8]      9.71
##    y0[9]     10.98
##    y0[10]     8.56
##    y0[11]     8.38
##    y0[12]    10.80
##    y0[13]    11.01
##    y0[14]     8.08
##    y0[15]     8.10
##    y0[16]     9.14
##    y0[17]    10.41
##    y0[18]     9.58
##    y0[19]     7.14
##    y0[20]    11.82
##    m1[1]      2.24
##    m1[2]      5.01
##    m1[3]      6.82
##    m1[4]      8.00
##    m1[5]      8.77
##    m1[6]      9.28
##    m1[7]      9.61
##    m1[8]      9.82
##    m1[9]      9.96
##    m1[10]    10.05
##    y1[1]      2.15
##    y1[2]      4.81
##    y1[3]      5.44
##    y1[4]      8.88
##    y1[5]      9.98
##    y1[6]      9.92
##    y1[7]     10.27
##    y1[8]      7.91
##    y1[9]     11.01
##    y1[10]    12.02
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -12.22 -11.89 1.34 1.06 -14.83 -10.80 1.01      632      847
##    b0      10.24  10.18 0.63 0.56   9.32  11.35 1.01      936      736
##    b1       0.54   0.53 0.12 0.10   0.36   0.75 1.01      742      624
##    s        1.27   1.24 0.23 0.21   0.95   1.70 1.00      845      929
##    m0[1]    6.01   6.01 0.56 0.53   5.10   6.95 1.01      784      708
##    m0[2]    7.11   7.12 0.51 0.49   6.24   7.93 1.01      825      768
##    m0[3]    9.84   9.85 0.40 0.39   9.19  10.48 1.00     1478     1335
##    m0[4]    7.10   7.10 0.52 0.49   6.23   7.92 1.01      824      768
##    m0[5]    3.41   3.38 0.47 0.43   2.70   4.23 1.01      755      622
##    m0[6]    9.01   9.01 0.33 0.32   8.46   9.54 1.00     1519     1161
##    m0[7]    2.31   2.28 0.35 0.32   1.79   2.94 1.01      750      640
##    m0[8]    9.41   9.42 0.32 0.32   8.87   9.93 1.00     2090     1357
##    m0[9]    9.84   9.85 0.40 0.39   9.19  10.49 1.00     1466     1279
##    m0[10]   9.17   9.17 0.32 0.31   8.63   9.69 1.00     1806     1422
##    m0[11]   8.77   8.78 0.35 0.34   8.19   9.32 1.00     1278     1009
##    m0[12]  10.02  10.02 0.47 0.45   9.27  10.83 1.00     1181     1081
##    m0[13]   9.94   9.95 0.44 0.42   9.25  10.67 1.00     1310     1213
##    m0[14]   9.95   9.96 0.44 0.42   9.25  10.69 1.00     1288     1149
##    m0[15]   9.26   9.27 0.32 0.31   8.73   9.78 1.00     2046     1468
##    m0[16]   8.72   8.73 0.35 0.35   8.13   9.27 1.00     1245      929
##    m0[17]   9.74   9.75 0.37 0.36   9.13  10.34 1.00     1649     1463
##    m0[18]   9.78   9.79 0.38 0.37   9.16  10.39 1.00     1586     1507
##    m0[19]   5.99   5.98 0.56 0.53   5.08   6.92 1.01      784      708
##    m0[20]   9.93   9.94 0.43 0.41   9.24  10.64 1.00     1331     1213
##    y0[1]    6.08   6.06 1.43 1.38   3.82   8.50 1.00     1847     1833
##    y0[2]    7.12   7.13 1.41 1.35   4.80   9.31 1.00     1693     1797
##    y0[3]    9.81   9.83 1.35 1.28   7.64  12.04 1.00     1879     1759
##    y0[4]    7.05   7.08 1.42 1.38   4.75   9.39 1.00     1584     1614
##    y0[5]    3.43   3.45 1.38 1.32   1.18   5.63 1.00     1879     1721
##    y0[6]    9.00   9.00 1.34 1.29   6.85  11.16 1.00     1846     1736
##    y0[7]    2.32   2.30 1.33 1.22   0.15   4.57 1.00     1804     1502
##    y0[8]    9.35   9.35 1.31 1.24   7.23  11.50 1.00     2082     1931
##    y0[9]    9.88   9.86 1.34 1.25   7.79  12.04 1.00     1828     1775
##    y0[10]   9.18   9.13 1.33 1.29   7.01  11.37 1.00     2305     1914
##    y0[11]   8.78   8.74 1.35 1.27   6.55  11.10 1.00     2016     1914
##    y0[12]  10.06  10.04 1.37 1.38   7.82  12.31 1.00     1628     1858
##    y0[13]   9.93   9.95 1.36 1.35   7.68  12.08 1.00     1883     1893
##    y0[14]   9.93   9.92 1.37 1.29   7.65  12.19 1.00     1958     1627
##    y0[15]   9.30   9.30 1.34 1.26   7.10  11.51 1.00     2156     1566
##    y0[16]   8.69   8.69 1.37 1.30   6.46  10.94 1.00     1837     1985
##    y0[17]   9.78   9.80 1.32 1.28   7.57  11.93 1.00     1731     1859
##    y0[18]   9.80   9.76 1.33 1.30   7.69  12.00 1.00     1956     2057
##    y0[19]   6.01   6.00 1.42 1.37   3.71   8.25 1.00     1554     1798
##    y0[20]   9.91   9.89 1.41 1.35   7.73  12.25 1.00     1849     1799
##    m1[1]    2.31   2.28 0.35 0.32   1.79   2.94 1.01      750      640
##    m1[2]    5.11   5.08 0.56 0.52   4.22   6.07 1.01      769      644
##    m1[3]    6.88   6.89 0.53 0.50   6.00   7.74 1.01      812      731
##    m1[4]    8.03   8.04 0.43 0.41   7.29   8.69 1.01      939      827
##    m1[5]    8.77   8.78 0.35 0.34   8.19   9.32 1.00     1278     1009
##    m1[6]    9.25   9.26 0.32 0.31   8.72   9.78 1.00     2038     1491
##    m1[7]    9.57   9.58 0.34 0.33   9.01  10.11 1.00     1933     1355
##    m1[8]    9.79   9.80 0.38 0.37   9.16  10.41 1.00     1567     1484
##    m1[9]    9.93   9.94 0.43 0.41   9.24  10.64 1.00     1329     1213
##    m1[10]  10.02  10.02 0.47 0.45   9.27  10.83 1.00     1181     1081
##    y1[1]    2.28   2.27 1.34 1.30   0.17   4.52 1.00     1892     1767
##    y1[2]    5.09   5.04 1.41 1.31   2.82   7.46 1.00     1610     1775
##    y1[3]    6.91   6.91 1.39 1.39   4.71   9.24 1.00     1617     1600
##    y1[4]    8.03   8.05 1.38 1.27   5.73  10.29 1.00     1791     1700
##    y1[5]    8.78   8.78 1.30 1.26   6.62  10.86 1.00     1953     1654
##    y1[6]    9.24   9.27 1.31 1.25   7.05  11.36 1.00     2051     1944
##    y1[7]    9.57   9.58 1.31 1.28   7.40  11.71 1.00     1876     1884
##    y1[8]    9.79   9.78 1.37 1.30   7.56  12.09 1.00     2093     1805
##    y1[9]    9.94   9.95 1.36 1.32   7.69  12.13 1.00     1700     1811
##    y1[10]  10.01   9.98 1.33 1.30   7.84  12.25 1.00     1664     1753

sigmoid curve

y=Ym/ 1+exp(-b1* (x-b0)) -> y~B(Ym, 1+exp(-b1*(x-b0)))
b0,b1>0
x[0,Infinity), y[0,Ym]
(x=b0, y=Ym/2)

ex8-5.stan

sigmoid curve

data{
  int N;
  vector[N] x;
  int Ym;
  array[N] int y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0,upper=100> b0;
  real<lower=0,upper=100> b1;
}
model{
  y~binomial_logit(Ym,b1*(x-b0));
}
generated quantities{
  array[N] int y0;
  for(i in 1:N){
    y0[i]=binomial_rng(Ym,inv_logit(b1*x[i]-b0*b1));
  }
  array[N1] int y1;
  for(i in 1:N1){
    y1[i]=binomial_rng(Ym,inv_logit(b1*x1[i]-b0*b1));
  }
}
n=20
x=runif(n,0,9)
ym=10
y=rbinom(n,ym,1/(1+exp(-2*(x-4))))
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,Ym=ym,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-5.stan')

mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "y0"   "y1"
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -44.74 -44.44 1.01 0.71 -46.81 -43.80 1.00      456      914
##    b0       3.86   3.86 0.12 0.12   3.67   4.07 1.00     1195     1004
##    b1       2.39   2.33 0.49 0.45   1.70   3.33 1.01      352      340
##    y0[1]    9.96  10.00 0.21 0.00  10.00  10.00 1.00     1655       NA
##    y0[2]    4.46   4.00 1.70 1.48   2.00   7.00 1.00     2055     1935
##    y0[3]    9.67  10.00 0.63 0.00   8.00  10.00 1.00     1240       NA
##    y0[4]   10.00  10.00 0.00 0.00  10.00  10.00   NA       NA       NA
##    y0[5]   10.00  10.00 0.05 0.00  10.00  10.00 1.00     2027       NA
##    y0[6]    5.38   5.00 1.74 1.48   3.00   8.00 1.01     1430     1911
##    y0[7]    7.40   8.00 1.58 1.48   5.00  10.00 1.00     1337       NA
##    y0[8]    3.05   3.00 1.54 1.48   1.00   6.00 1.00     1842     1860
##    y0[9]    2.31   2.00 1.41 1.48   0.00   5.00 1.00     1853     1902
##    y0[10]   0.01   0.00 0.12 0.00   0.00   0.00 1.00     1834     1834
##    y0[11]   0.44   0.00 0.66 0.00   0.00   2.00 1.00     1881     1915
##    y0[12]  10.00  10.00 0.06 0.00  10.00  10.00 1.00     2025       NA
##    y0[13]  10.00  10.00 0.03 0.00  10.00  10.00 1.00     2018       NA
##    y0[14]  10.00  10.00 0.04 0.00  10.00  10.00 1.00     2019       NA
##    y0[15]   1.74   2.00 1.28 1.48   0.00   4.00 1.00     1505     1841
##    y0[16]   0.04   0.00 0.21 0.00   0.00   0.00 1.00     1803     1774
##    y0[17]   0.01   0.00 0.10 0.00   0.00   0.00 1.00     2057     2057
##    y0[18]   2.28   2.00 1.42 1.48   0.00   5.00 1.00     2010     1970
##    y0[19]   0.31   0.00 0.58 0.00   0.00   1.00 1.00     1712     1693
##    y0[20]   0.03   0.00 0.17 0.00   0.00   0.00 1.00     1931     1954
##    y1[1]    0.01   0.00 0.11 0.00   0.00   0.00 1.00     2065     2065
##    y1[2]    0.07   0.00 0.27 0.00   0.00   1.00 1.00     1954     1972
##    y1[3]    0.41   0.00 0.66 0.00   0.00   2.00 1.00     1389     1413
##    y1[4]    2.58   2.00 1.46 1.48   0.00   5.00 1.00     1739     1962
##    y1[5]    7.39   8.00 1.57 1.48   5.00  10.00 1.00     1404       NA
##    y1[6]    9.52  10.00 0.77 0.00   8.00  10.00 1.00     1118       NA
##    y1[7]    9.92  10.00 0.30 0.00   9.00  10.00 1.00     1748       NA
##    y1[8]    9.98  10.00 0.12 0.00  10.00  10.00 1.00     2080       NA
##    y1[9]    9.99  10.00 0.07 0.00  10.00  10.00 1.00     1267       NA
##    y1[10]  10.00  10.00 0.04 0.00  10.00  10.00 1.00     2020       NA

y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,ymed=smy1$median,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=ymed),xy,col='black')

up and down

y=b0(exp(-b1* x)-exp(-b2* x)) -> y~N(b0(exp(-b1* x)-exp(-b2* x)),s)
b0,b1,b2>0, b1<b2
x[0,Infinity), 0<y<b0
(x=log b1-log b2 / b1-b2, y=max(y))

ex8-6.stan

up and down

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real<lower=0,upper=200> b0;
  real<lower=0,upper=1> b1;
  real<lower=0,upper=1> b2;
  real<lower=0,upper=10> s;
}
model{
  y~normal(b0*(exp(-b1*x)-exp(-b2*x)),s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0*(exp(-b1*x[i])-exp(-b2*x[i]));
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0*(exp(-b1*x1[i])-exp(-b2*x1[i]));
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,50)
y=rnorm(n,100*(exp(-0.03*x)-exp(-0.2*x)),1)
qplot(x,y)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-6.stan')

fn(mdl,data)
## Initial log joint probability = -5547.9 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
##       54      -7.43127   1.90425e-05   0.000920521           1           1      118    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__      -7.43
##    b0        95.89
##    b1         0.03
##    b2         0.21
##    s          0.88
##    m0[1]     56.55
##    m0[2]     51.35
##    m0[3]     45.09
##    m0[4]     25.45
##    m0[5]     60.17
##    m0[6]     48.96
##    m0[7]     30.89
##    m0[8]     30.35
##    m0[9]     60.43
##    m0[10]    54.36
##    m0[11]    30.60
##    m0[12]    58.95
##    m0[13]    59.20
##    m0[14]    53.62
##    m0[15]    40.88
##    m0[16]    36.12
##    m0[17]    60.56
##    m0[18]    27.90
##    m0[19]    49.48
##    m0[20]    26.63
##    y0[1]     54.94
##    y0[2]     50.94
##    y0[3]     46.29
##    y0[4]     26.44
##    y0[5]     59.58
##    y0[6]     48.06
##    y0[7]     30.93
##    y0[8]     30.34
##    y0[9]     60.29
##    y0[10]    53.52
##    y0[11]    31.04
##    y0[12]    58.90
##    y0[13]    58.66
##    y0[14]    53.37
##    y0[15]    43.07
##    y0[16]    34.99
##    y0[17]    60.97
##    y0[18]    26.25
##    y0[19]    50.06
##    y0[20]    26.61
##    m1[1]     54.36
##    m1[2]     60.57
##    m1[3]     58.27
##    m1[4]     53.23
##    m1[5]     47.60
##    m1[6]     42.19
##    m1[7]     37.25
##    m1[8]     32.83
##    m1[9]     28.91
##    m1[10]    25.45
##    y1[1]     54.67
##    y1[2]     60.44
##    y1[3]     59.99
##    y1[4]     53.23
##    y1[5]     46.92
##    y1[6]     42.42
##    y1[7]     36.12
##    y1[8]     32.95
##    y1[9]     28.47
##    y1[10]    23.91
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.2 seconds.
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.3 seconds.
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -11.31 -11.03 1.43 1.36 -13.98  -9.54 1.01      472      813
##    b0      95.93  95.80 2.76 2.53  91.68 100.52 1.01      495      738
##    b1       0.03   0.03 0.00 0.00   0.03   0.03 1.01      588      838
##    b2       0.21   0.21 0.01 0.01   0.20   0.23 1.01      525      773
##    s        1.01   1.00 0.19 0.18   0.70   1.35 1.15       21       28
##    m0[1]   56.54  56.53 0.37 0.37  55.92  57.13 1.02      720      971
##    m0[2]   51.34  51.35 0.36 0.35  50.76  51.93 1.01      680     1006
##    m0[3]   45.08  45.09 0.30 0.30  44.58  45.57 1.01     1546     1293
##    m0[4]   25.46  25.46 0.41 0.40  24.78  26.14 1.01      897     1365
##    m0[5]   60.15  60.15 0.43 0.40  59.45  60.88 1.02     1706     1646
##    m0[6]   48.96  48.96 0.34 0.33  48.40  49.51 1.01      816     1185
##    m0[7]   30.90  30.90 0.34 0.33  30.33  31.46 1.01     1291     1564
##    m0[8]   30.36  30.36 0.35 0.34  29.78  30.93 1.01     1235     1600
##    m0[9]   60.41  60.41 0.35 0.34  59.84  60.98 1.02     1885     1322
##    m0[10]  54.37  54.36 0.71 0.67  53.20  55.54 1.01      842     1342
##    m0[11]  30.61  30.62 0.35 0.34  30.04  31.18 1.01     1260     1600
##    m0[12]  58.93  58.92 0.36 0.36  58.34  59.51 1.02     1114     1045
##    m0[13]  59.20  59.19 0.50 0.49  58.36  60.04 1.02     1443     1386
##    m0[14]  53.60  53.60 0.37 0.37  52.98  54.21 1.02      659     1054
##    m0[15]  40.88  40.87 0.28 0.27  40.42  41.33 1.01     1997     1639
##    m0[16]  36.12  36.13 0.29 0.29  35.65  36.60 1.01     1833     1846
##    m0[17]  60.55  60.54 0.36 0.35  59.96  61.13 1.01     1972     1395
##    m0[18]  27.91  27.91 0.38 0.36  27.28  28.53 1.01      987     1576
##    m0[19]  49.47  49.47 0.34 0.34  48.91  50.02 1.01      771     1162
##    m0[20]  26.64  26.64 0.40 0.38  25.98  27.29 1.01      938     1368
##    y0[1]   56.52  56.53 1.13 1.09  54.68  58.30 1.01     1868     1718
##    y0[2]   51.32  51.32 1.11 1.08  49.50  53.09 1.00     1762     1874
##    y0[3]   45.08  45.08 1.06 0.99  43.32  46.84 1.01     1934     1805
##    y0[4]   25.43  25.43 1.12 1.06  23.59  27.22 1.00     1676     1370
##    y0[5]   60.12  60.15 1.09 1.04  58.29  61.86 1.00     1890     1798
##    y0[6]   48.95  48.97 1.14 1.08  47.05  50.79 1.00     2008     1749
##    y0[7]   30.89  30.86 1.09 1.04  29.11  32.70 1.01     2033     1731
##    y0[8]   30.36  30.39 1.09 1.07  28.56  32.08 1.01     1889     1775
##    y0[9]   60.44  60.43 1.06 1.00  58.78  62.17 1.01     1598     1619
##    y0[10]  54.35  54.31 1.25 1.20  52.29  56.40 1.00     1514     1543
##    y0[11]  30.63  30.63 1.13 1.10  28.84  32.46 1.00     1828     1769
##    y0[12]  58.97  58.96 1.11 1.11  57.14  60.78 1.01     2078     1585
##    y0[13]  59.18  59.19 1.15 1.07  57.21  61.11 1.01     1768     1663
##    y0[14]  53.65  53.64 1.12 1.05  51.80  55.54 1.00     2040     1696
##    y0[15]  40.87  40.90 1.06 1.03  39.11  42.59 1.00     1693     1585
##    y0[16]  36.12  36.14 1.07 1.04  34.38  37.84 1.01     1848     1662
##    y0[17]  60.57  60.55 1.08 1.07  58.88  62.40 1.01     1708     1865
##    y0[18]  27.91  27.89 1.09 1.06  26.13  29.69 1.00     1801     1808
##    y0[19]  49.46  49.44 1.08 1.05  47.68  51.35 1.01     1932     1823
##    y0[20]  26.65  26.67 1.10 1.04  24.81  28.40 1.00     1896     1677
##    m1[1]   54.37  54.36 0.71 0.67  53.20  55.54 1.01      842     1342
##    m1[2]   60.56  60.56 0.38 0.36  59.93  61.18 1.01     1960     1533
##    m1[3]   58.26  58.24 0.36 0.36  57.66  58.84 1.02      929     1004
##    m1[4]   53.21  53.21 0.37 0.36  52.60  53.82 1.02      658      923
##    m1[5]   47.59  47.59 0.33 0.32  47.06  48.13 1.01      963     1211
##    m1[6]   42.19  42.19 0.28 0.28  41.71  42.65 1.01     1882     1497
##    m1[7]   37.25  37.25 0.28 0.28  36.78  37.71 1.02     1942     1809
##    m1[8]   32.83  32.83 0.32 0.31  32.30  33.37 1.01     1472     1658
##    m1[9]   28.92  28.92 0.37 0.36  28.31  29.52 1.01     1086     1571
##    m1[10]  25.46  25.46 0.41 0.40  24.78  26.14 1.01      897     1365
##    y1[1]   54.37  54.37 1.20 1.16  52.43  56.30 1.01     1676     1704
##    y1[2]   60.58  60.58 1.09 1.04  58.82  62.36 1.01     2040     1825
##    y1[3]   58.26  58.26 1.10 1.07  56.47  60.05 1.01     1955     1357
##    y1[4]   53.21  53.20 1.05 1.01  51.53  54.94 1.01     2089     1792
##    y1[5]   47.60  47.59 1.05 1.00  45.87  49.29 1.00     1895     1683
##    y1[6]   42.17  42.16 1.06 1.00  40.50  43.90 1.01     2036     1750
##    y1[7]   37.26  37.26 1.12 1.08  35.40  39.10 1.00     1977     1361
##    y1[8]   32.80  32.81 1.08 1.02  31.03  34.51 1.00     1928     1768
##    y1[9]   28.91  28.92 1.12 1.09  27.03  30.64 1.01     1773     1678
##    y1[10]  25.46  25.46 1.09 1.03  23.61  27.29 1.01     1537     1648